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Abstract 


Solutions for transonic viscous and inviscid flows using a composite 
velocity procedure are presented. The velocity components of the 
compressible flow equations are written in terms of a multiplicative 
composite consisting of a viscous or rotational velocity and an 
inviscid, irrotational, potential-like function. This provides for an 
efficient solution procedure that is locally representative of both 
asymptotic inviscid and boundary layer theories. A modified 
conservative form of the axial momentum equation that is required to 
obtain rotational solutions in the inviscid region is presented and a 
combined conservation/non-conservation form is applied for evaluation of 
the reduced Navier-Stokes(RNS) , Euler and potential equations. A 
variety of results are presented and the effects of the approximations 
on entropy production, shock capturing and viscous interaction are 
discussed. 

1 . Introduction 

The composite velocity formulation developed by Rubin and Khoslatl] 
is a boundary layer like relaxation procedure based on a multiplicative 
composite velocity. For the Navi er-Stokes or reduced Navier-Stokes 
(RNS) equations it provides a technique that is consistent with both 
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asymptotic inviscid theory and boundary layer theory. The composite 
velocity formulation is directly applicable to incompressible, subsonic, 
transonic and supersonic flows. The solution procedure is identical in 
all cases. Supersonic regions are modelled with the Enquist*^Osher 
approximation. 

In this formulation a multiplicative composite for the velocities is 
defined, see reference (1). A viscous or rotational velocity (U) and a 
pseudo-potential or irrotational velocity (<|> x , 4yK are specified. This 

allows for the application of procedures developed for the transonic 
potential equation to be adapted to a composite velocity Euler and RNS 
formulation. Finally, the composite splitting provides for greater 
flexibility, since the viscous and inviscid portions of the velocity 
fields are easily Identified. 

The composite velocity formulation has previously been used to 
determine incompressible [2], subsonic and transonic flows for 
boattail geometries. Both laminar and turbulent flows were considered. 
In these studies a coupled (U-$) strongly implicit procedure [5] was 
used to solve the full Navier-* Stokes equations. In the present study 
application of the composite velocity technique for both viscous and 
inviscid transonic flows are presented. For viscous flows a reduced 
form of the compressible Navier-Stokes equations, in which the viscous 
terms in the normal momentum equation and the streamwise diffusion terms 
in the axial momentum equation are neglected, is considered. The 
continuity and streamwise momentum equations are solved using a coupled 
line relaxation procedure allowing for interaction between the boundary 
layer and the outer inviscid flow. The Enquist-Osher flux biasing 
scheme for the transonic potential equation is incorporated into the 
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solution procedure. Finally, a Cebeci**Smith two layer model [6] is used 
for turbulence closure. This model is acceptable for the attached or 
mildly separated wall layers considered herein. 

Verification of solution accuracy, in the outer inviscid region of a 
viscous interacting flow field, has been considered by solving the 
inviscid Euler equations. The composite velocity equations are still 
specified, but with slip boundary conditions and infinite Reynolds 
number (Re). Both irrotational, isentropic (potential) solutions and 
rotational, nonisentropic (Euler) solutions may be calculated with the 
Re - • and slip boundary conditions. A conservation form of the 
composite velocity equations is required for transonic Euler solutions; 
a non-conservation form generates the full potential solution. The 
conservation form for the Euler equations is shown to generate a correct 
entropy rise at the shock wave, while the non-conservation form does not 
generate spurious entropy in non shock regions. 

2. Governing Equations 

The Navier-*Stokes equations for steady, compressible flow of a 
perfect gas are given in general, orthogonal, curvilinear coordinates as 
follows: 

Continuity 

(ph 2 h^u)^ + (ph.jhgV)^ - 0 (U 

E-Momentum 


Z puu ♦ ^ pvu ♦ —-r-puvh, - r~— pv*h - - r P r + 

V 5 h 2 n h 1 h 2 1 n h 1 2 2 S h < 5 


( 2 ) 
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n-Momenturo 
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(3) 


DUc (h 2 h 3 T t2 } * 3n (h 1 h 3 T 22^ * T 1 


♦ T 


2 h t h 2 ' *11 * T 33 h 2 h 3 


Energy 


F, puH c * S » ,H n • ‘ 5F (h 3 h t q 2 )] * * ' T 
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where ♦ and f are the dissipation and viscous work terms , respectively 
and 


T 11 “ 2u ^h 1 u C * hjh^lJ "* 3 D^ (h 2 h 3 u) e + (h 3 h l v) n^ 
t 22 " 2y ^h7 V C * hlK?! J * 3 D^ (h 2 h 3 u) ^ + (h 3 h 1 v, J 

I \ C. ^ 

T 33 " 2w ^h 3 h 1 h 3 n + V^ h 3 n ^ "* 3 D^ (h 2 h 3 u) t + (h 3 h 1 v) J 


t 21 ” t 12 " ^h^h^C + h 2 ( h 1 ) n^ 


D - h 1 h 2 h 3 


1 . 3T 1 . 3T . _ 

q 1 ‘ h7 k 85 ! q 2 ’ ~ k 5F ! »• k * T 


In these equations, £ and u are the coordinate and velocity measured 
along the body surface; n and v are the coordinate and velocity normal 
to the body surface; p is the density, p the pressure, T the 
temperature, and H the total enthalpy. The terms h^(£,n), h 2 (£,n), and 


h 3 are the metrics for the curvilinear coordinate system. Two 

additional state equations are required to complete the governing set of 
equations. These are 
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where S is entropy. 

The composite velocity formulation procedure developed by Rubin and 
Khosla[1] is employed to represent the velocity fields. In the spirit 
of matched asymptotic expansions, the velocities are re-written as 

u - ""- - - (Hi) - (U+1 )u (7a) 

n 1 t e 


The multiplicative composite that represents the axial velocity consists 
of two terms, an irrotational "pseudo" potential function and a viscous 
velocity U. Since the change in v across the boundary layer is of the 
order of the boundary layer thickness, the normal velocity is determined 
solely by the "pseudo" potential function. 

By substituting equations (7 a,b) into the Navier-Stokes equations 
(1-3), the following system for U, $, S for 2-D (h^-l) conformal (i^hj) 

coordinates is obtained: 

Continuity 

[p(U+1 )(1+4>^)] 5 ♦ [P4> n ] n “ 0 
^-Momentum 


i[(ph 2 (U 2 *U)u e *)|r + ( P h i Uu e v) J + h^Ve^ 

-pu vUh. - ~S r - § G. ♦ viscous terms 
D k e 1 h t C h 1 ^ 
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n-Momentum 





The new variable G, that appears in these equations Is similar to the 
total (or Bernoulli-like) pressure. G is not, however, assumed to be 
constant, but is calculated by the solution procedure. The entropy S is 
zero in the undisturbed flow. 


In the viscous region, the composite velocity formulation is 
representati ve of the full or reduced form of the Navier M Stokes 
equations. The continuity and (^momentum equations determine 4 and U 
and the viscous total pressure correction is determined from the n~ 
momentum equation. In the limit as U+0 the continuity equation reduces 
to the full potential equation and the Bernoulli relation, G-constant, 
is recovered. Thus, equations (9) and (10) are identically satisfied 
and the composite velocity system has reduced to the expected 
representation .for an inviscid, irrotational flow. 

Inviscid flows are solved by dropping the viscous terms in equations 
(9) and (10). The interpretation of the composite velocity terms for 
inviscid flows varies slightly from that for viscous flows. The 4 term 
still represents an irrotational "pseudo" potential function. The U 
term, however, is no longer associated with the viscous effects but 
rather it reflects with the effects of rotat tonality associated with the 
inviscid flow; viz, the vorticity 0 - ^ h 2 V ^~^1 U ^n " “^(l-*-^))^. 
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The composite velocity scheme is formulated to provide solutions to 
the full Euler equations; if the system is solved in the non-* 
conservation form, however, the full potential solution is recovered 
instead; i.e., U-0 everywhere and neither entropy nor vorticity is 
generated at any point in the flow field. This includes leading and 
trailing edges as well as captured shock waves. In order to capture the 
rotational or (Euler) shock wave in transonic flows, the ^-momentum 
equation must be rewritten in a quasi-conservation form. This is 
obtained by substituting for S and G in the right hand side of equation 

u 

(9), i.e., subtract, from the right hand side, r— r~ times the continuity 

n l 2 

equation. This gives the following ’conservation' form for the 
momentum equation: 

jj[(ph 2 (U 2 +2U)u*) e + (ph 1 Uu e v) n J ♦ ^PU e vUh 1 - (11) 


h 2 h 1 

-£-<r *pu>) e - 1- (p V ) n * ^(pv* " pu’> - 2.o-53<p V ). 

For subsonic flows (U*0) and a cartesian grid (hj-h 2 *1.0), equation 


(11) reduces to the familiar conservation form of the ^-momentum 
equation, i.e., 

(P ♦ pu 2 )_ + (pu v) - 0 (12) 

e £ e n 

For transonic flow the correct entropy rise at the shock wave will now 
be generated. Although the correct entropy rise is predicted at the 
shock with the system (11), spurious entropy is also generated in non- 
shock regions. The generation of numerical entropy is a common problem 
found in many Euler solvers. Large errors in entropy may be generated 
at leading and trailing edges[7]. These errors may even lead to 


7 


spurious unsteady or steady solutions [8]. In the present technique a 
simple solution to this problem is available. The nonconservative form 
of the axial momentum equation (9) produces no entropy, but will 
accurately convect entropy or vorticity that is generated elsewhere. 
Therefore, this form of the axial momentum equation is used everywhere 
except in the shock region, where equation (11) is required. This leads 
to a solution procedure with the desirable feature of generating the 
correct entropy rise at the shock wave, but not creating spurious 
entropy in other regions of the flow. The advantage of non M conservative 
equations away from shock waves and combined conservative/non- 
conservative systems has been discussed in several studies by other 
investigators[21 ]. 

3. Boundary Conditions 

As a simple test case, solutions for symmetric flow over a NACA0012 
airfoil are obtained. The boundary conditions for the composite 
velocity formulation are easily implemented. At the inflow, 

uniform flow is assumed; thus U»0, H*H , ♦=0, and S-0. The upper 

© 

boundary , n-n , is assumed to be sufficiently far from the airfoil so 
m 

that the flow field is undisturbed; therefore, similar boundary 
conditions apply. Along the body, the viscous no-*slip and zero 
injection conditions are used; therefore, U--1 and ♦ n *°* Ahead of the 

airfoil and in the wake, the symmetry conditions are imposed. 

At the outflow, £-£ , only ♦ -0 must be prescribed. This, in effect, 

m i, 

assumes a weak viscous/inviscid interaction. 

For inviscid flows, the no slip condition no longer applies and the 
zero vorticity condition U^-0 is specified along the airfoil. At the 
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outflow, £-£ m . the boundary condition ♦-O Is imposed. The remaining 

boundary conditions are unchanged. 

Finite Differencing 

The Enquist-Osher flux biasing scheme[9] for transonic flows has 
been adapted for differencing the continuity equation. This scheme, 
which has been developed for the full potential equation, produces very 
sharp shocks and guarantees that expansion shocks do not occur. This 

scheme consists of defining a modified density p, such that 


i.j ' p i,j "Sf.j <(p<,) - i.j ' (pq) ' i-i ,j> 


(13) 


where 

(pq)_ - 0.0 if M S 1 

* « 

(pq)_ - pq “ P q if M M 

x « 

Here p and q are the sonic velocity and density. The modified density 


p is then used in the [p(1 + ^)] portion of the ^-derivative in the 
continuity equation. 

The [p(U+1 )<!♦♦)] derivative is approximated with a two point 
backward difference and the derivative is approximated with a two 

point forward difference. The modified density p is applied as 
discussed previously. The C P<t> n ls differenced in a similar manner, 

when v is less than zero. When the v velocity is greater than zero, 
however, this term is forward differenced and is approximated with a 

two point backward difference. Central differencing results for 
subsonic regions. 
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In the ^-momentum equation (9), the (ph„(U*+U)u 2 ) derivative is 

approximated with two or three point backward differences. For viscous 
flows, upwind differencing is applied in regions of reversed flow. The 
convective term is Identically zero at the separation and reattachment 
points. The (ph^Uu^)^ and U g (or p^) terms are approximated with 

second order accurate central differences. For the and terms on 

the right hand side, two point backward differences are applied. When 

the conservation form, equation (11), is used, a two point backward 

difference is used for (p+pu 2 )_ and central differencing is used for 

® 5 

(pu v) . 
e n 

Backward differencing of the ^-derivatives provides for the proper 
convection of U in inviscid regions. If a value of U (i. e. vorticity) 
is generated, this value will be convected downstream in one global 
pass. Moreover, in the non-conservation form the value of the vorticity 
is conserved and additional numerical vorticity is not generated even at 
the trailing edge of an airfoil. 

The discretized continuity and 5-momentum equations are solved for U 
and <p with a coupled line relaxation procedure. The values of p, S, and 
G in the equations are given from the previous iteration. From the 
values of U and $, the entropy variation is determined from the n- 
momentum equation. This equation is solved using the standard 
trapezoidal rule. 

For the large Reynolds number moderately separated flows considered 
herein, a Cebeci-Smith two layer eddy viscosity model is used to close 
the system of governing equations. The coefficient of viscosity Is 
given as 
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where ^ la the laminar viscosity and ii fc la the turbulent eddy 

viscosity, as prescribed by the Cebeci-Smith model. Details of this 

closure model are found in reference 3. The onset of turbulence was 

♦ 

specified at 5K of the chord. The actual transition location is not 
well known, but is approximated reasonably with the present assumption. 

5. Results 

Viscous and inviscid solutions for transonic flow have been obtained 
for the symmetric flow over a NACA0012 airfoil. A Schwartz-Christoffel 
mapping procedure developed by Davis [10] generates the required grid, 
figure 1 . The flow region is defined from an upstream location of 
5— *1.6 to an outflow boundary 5-6.89. The airfoil is located between 
5-0.0 and 5-1.0. In the tangential directionj 105 grid points are 
prescribed; this includes 60 points on the body, 20 points ahead of the 
body, and 25 points aft of the body. A uniform mesh is defined on the 
body and a mesh stretch factor of 1.18 is used ahead of and aft of the 
body. In the n direction, there are 50 grid points, so that the flow 
region extends from n-0.0 to n-22.0. For viscous flows, an Initial grid 
spacing of An-0.0001 is prescribed; a stretch factor of 1.25 is assumed 
for the remaining mesh development. For inviscid flows, the initial 
grid spacing in the n-direction is An-0.02, since the boundary layer 
need not be resolved. 

To verify that the algorithm correctly calculates the inviscid 
portion of the flow, a full potential solution for H^-0.85 is obtained. 

Full potential solutions are obtained directly from the set of equations 
(8-10) by solving the axial momentum equation In the non-conservation 
form given by equation (9). In this form, U and S are calculated to be 


identically zero and therefore the full potential solution is recovered 
for The present solution is compared with results of the GAMM 
workshop on transonic flows [11], figure 2. The workshop results 
indicate a wide variation in shock location and strength and the figure 
depicts only the lower and upper bounds of the workshop solutions. 
Comparison with other potential solutions, in the published literature 
and for the same case, indicate that this scatter of solutions to the 
full potential equation is not unusual, figure 3> The present 
calculation produces a very sharp shock and the shock falls within the 
band of solutions presented at the GAMM workshop. 

Next, the alternate form of the axial momentum equation (12) is 
implemented in the shock region so that Euler solutions can be obtained. 
A comparison of a M^-0.8 solution, with results given by Clarke 

et.al.[19] is shown in figure 4. The solutions are seen to be in good 
agreement with the present results, which predict a slightly sharper 
shock. Figure 5 presents a comparison of the present solution with the 
GAMM workshop results for a M^-0.85 case. The band of solutions for the 

Euler case is seen to be much smaller than that for the full potential 
equation and again only the upper and lower bounds of solutions are 
presented. The present solution falls within this band and produces a 
sharper shock than those given in the workshop results. 

In figure 6, the entropy generated along the airfoil for the 
M *0.85 case is compared with the value obtained from the Rankine- 

Hugoniot shock relations. The agreement is seen to be excellent. Two 
other important features should be noted. First, entropy is not 
generated ahead of the shock region. Secondly, the entropy generated at 
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the shock is convected properly downstream with no additional increase 
in entropy elsewhere. 

The primary goal of this work, however, is to solve transonic, 
viscous interacting problems. The potential and Euler examples 
discussed herein have verified that the proper inviscid solutions are 
accurate, and therefore, should provide accurate outer flow behavior for 
full viscous calculations. The next set of results will describe 
solutions to the RNS equations for a variety of freestream Mach numbers 

and for a Reynolds number, Re-^xlO*. Comparison of these solutions for 
an irrotational (potential) outer inviscid flow model and for a 
rotational (Euler) outer model flow is included in the discussion. 

Results for M -0.8 are presented in figures 7a-c. In figure 7a the 

pressure coefficient is shown. Little difference is seen between the 
solutions with potential or Euler outer flows models. These results are 
also compared with experimental data [20]. The computed solution shows 
good agreement with the experimental results, except that the shock lies 
slightly aft of the experimental data. The skin friction coefficient 
for both outer flow conditions is compared in figure 7b. The effect of 
the outer flow is much more pronounced for the skin friction than was 
the case for the pressure coefficient. The skin friction has a much 
smaller decrease through the shock for the Euler outer flow. The 
oscillations in the skin friction coefficient at the leading edge are 
due to a lack of grid resolution and vanish at 10% of the chord. 
Finally, the Mach contours for M^-0.8 are given in figure 7c. 

Results for M -0.83 are presented in figures 8 a-c. The pressure 

oo 

coefficient, figure 8a, is once again insensitive to the potential or 
Euler outer flow modeling, except in the post shock region where the 
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Euler outer flow provides a somewhat different character than that 
obtained with the potential outer flow. The computed results again show 
the shock located slightly aft of the experimental data. Also, the 
experimental results 3how a small increase in pressure ahead of the 
shock, which is not seen in the computed results. Only the Euler outer 
flow model provides the proper flow character after the shock. The 
reason for the difference in solutions after the shock can be seen from 
the skin friction results, figure 8b. The solution with a potential 
outer flow model has a small separated region after the shock. The 
solution with an Euler outer flow model does not separate and the skin 
friction is somewhat larger. Mach contours for the M^-0.83 case are 

given in figure 8c. 

The importance of the outer flow model becomes more apparent for the 
Mach number M =0.85. These results are given in figures 9a-c. The 

shock obtained with the potential outer flow model lies ahead of the 
shock obtained with the Euler outer flow model and again predicts a 
different post shock behavior. The skin friction coefficient figure 9b, 
shows that the Euler outer flow model leads to a rapid recovery after 
the shock with only a small separation region; the potential outer flow 
model produces a large region of separation. The rotational ity in the 
outer flow appears to suppress the tendency toward separation. The Mach 
contour plot is given in figure 9c. Experimental results are not 
available for this case. 

The potential, Euler, and RNS (Euler outer model) results for M^-0.8 
and M^-0.85 are compared in figures 10 and 11, respectively. For 
M -0.8, the three solutions agree fairly well, with the Euler and RNS 

09 
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solution being slightly weaker and lying forward of the potential shock. 
The rotational and viscous effects are of increased importance for the 
larger Mach number M^-0.85. The Euler shock is weaker and lies forward 

of the potential shock. The RNS shock lies forward and is slightly 
weaker than the Euler shock. For the larger Mach numbers, shocks 
obtained from the full potential solution are too strong and are located 
far downstream of the RNS (Euler outer flow) results. 

6. Summary 

The composite velocity solution procedure has been applied for both 
viscous and inviscid transonic flows. Results are presented for flow 
over a NACA0012 airfoil. The results demonstrate the versatility of the 
composite velocity procedure. Both viscous and inviscid flows are 
solved from the same formulation with a simple change in the boundary 
conditions. 

For inviscid flows, both irrotational potential solutions and 
rotational Euler solutions are obtained. For the Euler model, the axial 
momentum equation is given in a full conservation form in the shock 
region. This provides a solution technique that produces the correct 
entropy rise at the shock and at the same time convects the entropy 
accurately; no spurious entropy is created outside of shock regions. The 
potential and Euler solutions were found to agree with earlier results 
presented for the inviscid models. 

Solutions for the RNS equations are obtained with potential and 
Euler outer models for a variety of Mach numbers. The results agree 
quite well with experimental results. The form of the outer flow model 
affects the post shock solution. Solutions with an irrotational outer 
flow tend to more readily induce post-shock separation, than do the 
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solutions with rotational outer flow modelling. The ability of the 
composite velocity procedure to efficiently calculate high Reynolds 
number transonic viscous flows attributes to the robustness of the 
solution technique. 

Finally, the results for potential, Euler, and RNS (Euler outer 
modelling) solutions for several cases are compared. For higher Mach 
numbers, the Euler and potential solutions produce stronger shocks that 
are located further aft on the airfoil and do not accurately reflect the 
shock behavior. 
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Figure 1. Viscous Grid for on NACA0Q12 Airfoil 
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Figure 2. Comparison of Potential Solution for an 
NACAO0 1 2 Airfoil, M — O.SS 



Figure 3. Comparison of Present Potential 

Solution with Various other Potential 
Solutions for an NACA0012 Airfoil 



Figure 4. Comparison of Euler Solution 

for an NACA0012 Airfoil. M = 0.8 





Figure 6. Comparison of Entropy along an 
NACA0012 Airfoil with the 
Rankine— Hugoniot Value, K/1 *= 0.85 
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with Different Outer Flows 



JO 


with Different Outer Flows 





Figure 7c. Mach Contours for an NACA0012 

Airfoil, M^-O.8, Re=*4.0x10 e ; 

R.N.S. Solution with Euler Outer Flow 
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with Different Outer Flows 
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Figure 8b. R.N.S. Solution for an NACA0012 

Airfoil, M oo =0.83, Re=4.0x10 e ; 
Comparison with Different Outer Flows 



Figure 8c. Mach Contours for an NACA0012 

Airfoil. M oo =»0.83, Re»4.Ox1O 0 ; 

R.N.S. Solution with Euler Outer Flow 



Figure 9a. R.N.S. Solution for an NACA0012 

Airfoil, M^O.85, Re=4.0x10 a ; 
Comparison with Different Outer Flows 



Figure 9c. Mach Contours for an NACA0012 

Airfoil, M^asO.85, Re=»4.0x10*; 

R.N.S. Solution with Euler Outer Flow 








